Vortex solitons in dipolar Bose-Einstein Condensates 
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We predict solitary vortices in quasi-planar condensates of dipolar atoms, polarized parallel to 
the confinement direction, with the effective sign of the dipole-dipole interaction inverted by means 
of a rapidly rotating field. Energy minima corresponding to vortex solitons with topological charges 
i = 1 and 2 are predicted for moderately strong dipole-dipole interaction, using an axisymmetric 
Gaussian ansatz. The stability of the solitons with ^ = 1 is confirmed by full 3D simulations, whereas 
their counterparts with i = 2 are found to be unstable against splitting into a set of four fragments 
(quadrupole). 

PACS numbers: 



I. INTRODUCTION 

Matter-wave patterns in Bose-Einstein condensates 
(BECs) are sustained by the interplay between the exter- 
nal trapping potential and intrinsic interactions between 
atoms. In particular, bright solitons and soliton trains in 
nearly one-dimensional (ID) traps [1] are supported by 
the relatively weak attraction between atoms of ^Li, or 
the stronger attraction in the ^^Rb condensate. For re- 
pulsive interactions, adding an axial optical-lattice (OL) 
potential gives rise to gap solitons, as demonstrated in 

In planar 2D geometry with intrinsic repulsion, delo- 
calized vortices constitute basic BEC patterns [3]. The 
creation of 2D matter-wave solitons (as well as 2D spa- 
tiotemporal solitons, alias "light bullets", in nonlinear 
optics [4]) is a challenge, as the contact attraction leads 
to collapse in this case. Square-shaped OL can stabi- 
lize fundamental solitons and solitary vortices (solitons 
with embedded vorticity) in two dimensions 0, @. As 
concerns vortex solitons, the lattice breaks the isotropy 
of the embedding space and the related angular momen- 
tum conservation. Nevertheless, the intrinsic topological 
charge of the vortex can be defined unambiguously as 
£ = A<I>/27r, where A<I> is the circulation of the phase of 
the respective complex wave function around the vortex' 
pivot. The simplest "crater-shaped" vortices, in the form 
of a single density peak with an inner hole, trapped, es- 
sentially, in a single cell of the square lattice, is unstable 
0. Stable vortices with charge £ can be constructed, in 
the simplest form, as sets of four peaks, with phase shift 
Aip = 7r£/2 between adjacent ones. Two stable four- 
peak vortex structures are possible, "rhombuses" (with 
a nearly empty cell in the middle) [5] and more compact 
"squares" [6]. For £ = 2, i.e. A(p = tt, these patterns are 
actually quadrupoles. Higher-order stable vortices, up to 
^ = 6, were found too, in the form of circular chains of 
8 or 12 peaks [8]. Also found were "supervortices" , built 
as ring-shaped chains of 12 compact (crater-shaped) lo- 
cal vortices with individual vorticity Iq = 1^ onto which 
global vorticity ^ = ±1 is imprinted Q. The supervortex 
is stable despite the instability of crater-shaped vortices 



in isolation. Two-dimensional solitons, as well as vor- 
tices and quadrupoles of the rhombic type, can also be 
stabilized by the quasi- ID OLs, i.e., periodic potentials 
depending on a single coordinate [9]. 

For repulsive contact interactions, square-shaped OL 
can support 2D fundamental and vortical gap solitons 
[10]. In addition, axisymmetric radial potentials may 
stabilize solitons, including vortical ones, in both cases 
of the attractive and repulsive interactions [llj . 

Despite the theoretical progress, 2D matter-wave soli- 
tons have not yet been observed, vortex solitons being a 
still more challenging subject. Therefore, the search for 
viable settings allowing the realization of such structures 
remains highly relevant. In parallel to BEC, theoretical 
and experimental studies of multidimensional spatiotem- 
poral solitons draw great interest in nonlinear optics Q. 

New possibilities for producing 2D solitons emerge in 
dipolar quantum gases, such as BECs of magnetically po- 
larized ^^Cr atoms [l2j, dipolar molecules [13], or atoms 
with electric moments induced by strong electric field 
[13] or laser illumination [l5[. Of these systems, the gas 
of chromium atoms is the medium available for current 
experiments. For axisymmetric geometry, with dipoles 
polarized perpendicular to the 2D plane, the natural 
anisotropic dipole-dipole (DD) interaction gives rise to 
in-plane repulsion and axial attraction, which can sup- 
port vortex lattices [16] and, in principle, 2D gap soli- 
tons (in the presence of the respective OL). On the other 
hand, the sign of the DD interaction may be reversed 
by means of rapid rotation of the dipoles \17v or using a 
combination of microwave and dc fields^, which en- 
ables the creation of isotropic solitons [ly]. However, in 
the full 3D geometry, isotropic vortex lines are destabi- 
lized by the DD interactions [20]. Alternatively, stable 
anisotropic solitons can be supported by the natural DD 
interaction, when dipoles are polarized in the 2D plane 
[21]. Related work in nonlinear optics employed the non- 
local thermal nonlinearity to predict stable vortex rings, 
with topological charges £ = 1 and 3 [22[ and 2D ellipti- 
cally shaped spatial solitons [2O]. 

Here, we assume the axisymmetric configuration as in 
Ref. [ly], with the aforementioned sign inversion of the 
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DD interaction fr?], to predict stable 2D vortex solitons. 
Families of vortex soliton states, with i = 1 and 2, are 
constructed in the framework of the 3D Gross-Pitaevskii 
(GP) equation, including the long-range DD interactions, 
contact repulsion, and transverse confinement potential. 
With topological charge ^ = 1, the solitary vortex is sta- 
ble, whereas the vortex soliton with i = 2 splits via a 
quadrupole instability. Since the only reported experi- 
mental demonstration of 2D spatial vortex solitons re- 
quires the presence of photonic lattices in photorefrac- 
tive crystals [25], and because no observation of vortex 
solitons was reported in uniform media, the proposal to 
create such solitons in dipolar BECs is pertinent to the 
experiment, especially in view of recent advances in tun- 
ing out the local nonlinearity via a Feshbach resonance 
(FR) [24]. 

The paper is organized as follows. In Section I, we 
apply a variational approximation to predict localized- 
vortex states yielding a minimum of the dipolar conden- 
sate energy. Naturally, these states have a chance to 
represent stable vortex solitons. In Section III, results 
of direct simulations of vortex solitons, performed in the 
framework of the full 3D GP equation, are summarized. 
For that purpose, we use both numerically exact pro- 
files, which are generated from the variational ansdtze 
by means of the preliminarily simulated propagation in 
imaginary time, and the ansdtze themselves, the corre- 
sponding results being quite similar. The paper is con- 
cluded by Section IV. In particular, in that section we 
discuss the physical significance of three-body losses in- 
duced by the FR. 



LJz is the transverse-trap frequency and Iz = ^/hjmuJz 
is its respective length. The interaction strengths are 

g = AirNas/lz and gd = NdPm (fi^lz) , where > is 
the 5-wave scattering length, d and m the atomic dipole 
moment and mass, N the number of atoms, and the nor- 
malization is taken in the form of j {y)\^ dr = 1. 

To approximate vortex-soliton states with topological 
charge we use the normalized Gaussian ansatz in cylin- 
drical coordinates p = Y^^^^"^r^, and (j) = tan~^(?//x), 

z, (j)) = Atp^ exp (- (ap^ + 72;^) /2) exp(i£^), (5) 

where A\ = {tt-^^^ /£)a^+^j^/'^ for ^ = 1,2. Evaluating 
El^pi}, we obtain. 
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where functions 
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/i(k) = -1 + 3/ R{K,x)[l + Q^{K,x)]dx, 
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h{K) = -1+3 
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II. VARIATIONAL ANALYSIS FOR ENERGY 
MINIMA 

We assume that a strong magnetic field aligns dipole 
moments along confinement axis z [13, [3, [3]. The 
respective GP energy functional, expressed in terms of 
BEG order-parameter (x, z), is 



(1) 



where the kinetic, confinement, and contact-interaction 
energies are, respectively. 



T=\j iV^WI^dr , V=\j z^\^{v)\Hv, 

U=\j |^(r)|^dr, 
and the DD mean-field energy is 
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Here and below, the length, time, and energy are scaled 
as r ^ y/Iz^ t uozt^ and E E/ {TiuOz)^ where 
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depend solely on the aspect ratio^ k, = \/^yJa. 

In strongly prolate (cigar-shaped) geometry with <C 
1, one has ^ 0, Q ^ 1, and /i,2(^) ^ — 1. By con- 
trast, for an oblate (pancake) shape with ^ 1, we have 

^ 1, Q ^ 0, and f 1,2(1^) 2. The change of the 
sign of /i,2 corresponding to the transition from <^ 1 
to /t: 1 is due to the respective change of the relative 
strength of "side-by-side" and "head-to-tail" DD interac- 
tions, which dominate the prolate and oblate configura- 
tions, respectively. Gonsequently, for gd > and fixed 7, 
the integrated DD energy, Ud, decreases both for a ^ 
(since Ud > for large hi) and a ^ 00 (because Ud < in 
this limit), leading to either expansion or collapse along 
P' 

Inversion of the sign of gd [HI, ll8|, ll9|, converts the 
energy maximum at the prolate/oblate transition point, 
into a minimum. It is thus required to stabilize 2D 
isotropic patterns, with the dipolar moments polarized 
along the cylindrical axis. This requirement for gd in- 
version, in combination with the necessity to reduce the 
strong contact repulsion using the FR, is an experimental 
challenge. The existence of stable anisotropic fundamen- 
tal (non-topological) solitons, with the dipoles polarized 



in the (x, y) plane and the ordinary sign of the DD inter- 
action {gd > 0) [21], suggests that anisotropic topological 
(vortex-like) solitons may be found in the same setting. 
However, we leave the analysis of such complex patterns 
to a separate work, aiming here to retain the cylindrical 
symmetry, adopting the assumption of < 0. 

Assuming the sign inversion of the DD interaction, the 
minimum in Ud as a function of a at 7 ~ 1 will trans- 
late into a minimum of total energy provided that 
Ud IS large enough to offset the contact-interaction and 
gradient terms, U and T, both scaling linearly with a. 
This requirement results in the necessary condition for 
the existence of a stable isotropic vortex soliton. 
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The inequalities on the right- and left-hand sides of Eq. 

guarantee, severally, dE{a^^) / da < for a ^ 
{fi 2) anda£;(a,7)/Sa > for a ^ oo (/^ ^ -1). As 
expected, conditions ([9|) can only be satisfied for gd < 0. 
The required strength, \gd\^ increases with i due to the 
growing centrifugal contribution to the kinetic energy, 
which must be balanced by the attractive part of the 
DD interaction. In the strong-interaction regime, with 
3^/47r ^ 6\^27r£, the kinetic term may be neglected, and 
Eq. ^ simplifies to -gd/g > 3/(87r). 

Condition (|9|) only guarantees a minimum of the GP 
energy as a function of a for fixed 7, which is not yet 
sufficient for a true minimum of E{a^ 7) in the (a, 7) 
plane. In particular, fixing and varying 7 = hz'^a^ one 
can see that U ^ Ud oc 7^/^, V ex I/7, and T ex 7, so 
that kinetic energy T cannot balance interaction terms 
U -\- Ud for large 7. Since U -\- Ud < in the IdLige-hz 
sector of the (a, 7) plane, this will lead to the 3D collapse, 
with a, 7 00, unless the vertical-confinement size Iz is 
smaller than the effective healing length, determined by 
the interplay of the kinetic energy with the combined 
contact and DD interactions. 

In Fig. 1, we plot the GP energies, as given by Eq. (|6]) 
(panels [T^-c) and Eq. ([7j) (panels [T]i-f), in three charac- 
teristic interaction-parameter regimes. For £ = 1 and 2 
respectively. Figs, la and Id display cases when the DD 
interaction is too weak to satisfy conditions (|9|). This re- 
sults in a monotonic decrease of the energy as a ^ at 
fixed 7, causing the radial expansion of the BEG. On the 
other hand, if the DD interaction is too strong, the energy 
decreases monotonically for fixed as a, 7 ^ 00, imply- 
ing the 3D collapse (Figs. Ic and If). In the intermedi- 
ate regime, the DD interaction is strong enough to offset 
the dispersive effect of the contact and kinetic terms, yet 
is not excessively strong to induce the 3D collapse. This 
regime gives rise to local energy minima, which are found 
at ^ = -gd = 20 for a = 0.081,7 = 1.15 in Fig. 1(b), 
and at ^ = 20, = -30 for a = 0.046,7 = 1.13 in Fig. 
1(e). These minima suggest the possibility of metastable 
oblate vortex solitons, with radial widths of 3.51/^ and 
4.66/2 for = 1,2 respectively. 



FIG. 1: (Color online) The GP energy functional for a vortex- 
soliton ansatz (|5j), with i = 1 (a-c) and i = 2 (d-f), as a 
function of variational parameters a and 7. In all plots, g = 
20, with gd = -10 (a), -20 (b), -30 (c) for i = 1, and 
= -20 (d), -30 (e), -40 (f), for i = 2. 



III. MEAN-FIELD DYNAMICS AND 
STABILITY 

To directly verify the existence and stability of axisym- 
metric vortex solitons, we substituted the values of a and 
7 corresponding to the local energy minima into ansatz 
([5]) , and simulated its evolution according to the full time- 
dependent GP equation. 
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First, the propagation in imaginary time was carried 
out to reshape the input into a numerically exact sta- 
tionary solitary vortex. The amplitude difference be- 
tween the initial Gaussian ansatz and reshaped soliton 
was < 10%. Then, to test the stability of the solitary 
vortices, we used these profiles as initial conditions and 
carried out 3D simulations in the real time. The lo- 
cal density and phase at 2; = 0, during the real-time 
evolution for approximately 20 trap periods, are shown 
in Figs. [2] and [3] for = 1, and in Figs. [3] and [5] for 
£ = 2. With topological charge ^ = 1, the vortex soli- 
ton remains virtually unchanged during the evolution 
(Fig. [2]), which demonstrates its full stability. By con- 
trast, with £ = 2^ the solitary vortex is modulationally 
unstable against splitting into a quadrupole set, as shown 
in Fig. m This observation is reminiscent of the stability 
analysis for solitary vortices in condensates with the lo- 
cal self-attraction, trapped in the axisymmetric parabolic 
potential where only £ = 1 vortices have their sta- 
bility region, all vortices with £ > 1 being inherently 
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FIG. 2: (Color online) The evolution of a vortex soliton with 
= 1, for g — —Qd — 20. The initial conditions are obtained 
by the imaginary-time propagation for icj^t = 30, starting 
from Gaussian ansatz (|5]) with a — 0.08 and 7 = 1.15, which 
corresponds to the energy minimum in Fig. lb. The evolution 
of the density profile in cross section z — demonstrates the 
stability of the solitary vortex with I — 1. 
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FIG. 4: (Color online) The evolution of the density profile 
in cross section z — for a vortex soliton with 1 — 2 and 
interaction strengths g — 20, Qd — —30. Initial conditions 
were generated by the imaginary- time propagation for iujzt — 
30, starting from Gaussian ansatz (|5]) with I — 2^ a — 0.046, 
7 = 1.13, which corresponds to the energy minimum in Fig. 
le. The solitary vortex is unstable, and eventually splits into 
a quadrupole set. 
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FIG. 3: (Color online) Phase plots during the evolution of the 
solitary vortex with i = 1 and the same parameters as in Fig. 

m 




FIG. 5: (Color online) Phase plots during the evolution of the 
solitary vortex with £ — 2 and the same parameters as in Fig. 

n 



unstable. 

The evolution of the solitary vortex with ^ = 1, was 
also compared to the CP dynamics in limit cases when 
either the contact repulsion or DD interaction is turned 
off. As mentioned above, the contact interaction may 
be experimentally controlled by means of the FR [24], 
and the DD interaction may be tuned using additional 
external fields [l7|, ll8|, or simply turned off by remov- 
ing the polarizing field, which allows the dipolar order 



to frustrate. In Fig. EJa-c), we display the density and 
phase distributions in cross section z = 0, for ^ = 1 and 
cozt = 12. Recall that the vortex soliton was robust for 
unaltered interaction strengths (Fig. [6^1). In contrast to 
that, when the DD interaction is turned off (Fig. [6)3), 
radial expansion is observed, being accompanied by an 
outgoing density current. On the other hand, when the 
gd/g ratio is too large, the vortex collapses (see Fig. [6t), 
and the current is funneled towards the vertical axis. It 
is noteworthy that the axial symmetry is maintained in 
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FIG. 6: (Color online) Density profiles (first line) and phase 
plots (second line) in cross section z = 0, after the real-time 
propagation for oozt = 12, starting from the vortex- soliton 
state with ^ = 1, as taken from Fig. 2. The robust evolution 
of the solitary vortex for g — —gd — 20 (a) is compared to 
the expansion for g = 2^^ ga — ^ (b) and collapse for ^ = 0, 
gd — —20 (c). Panels (d)-(f) show the real-time evolution of 
the density profile |'0(p, 2; = 0, 0)1^, starting with approximate 
Gaussian ansatz (O, for the same parameters. 



all cases. 

The stability of the solitary vortex with ^ = 1 is fur- 
ther corroborated by the simulated propagation in real 
time of the initial configuration corresponding to Gaus- 
sian ansatz ([5]), with a and 7 set to the location of the 
minimum in Fig. lb. This state, which actually cor- 
responds to a slightly perturbed vortex soliton, demon- 
strates, for g = —gd = 20, stable self-trapping into the 
exact solitary vortex (Fig. [6li). By contrast, expansion 
(Fig. [6^) and collapse (Fig.[6f) are observed, respectively. 



in the regimes of too weak (as in Fig. [T^) and too strong 
(Fig. [It) DD interactions. 



IV. CONCLUSIONS 

The experimental feasibility of quasi-2D soHtons in 
dipolar BECs was estimated in Refs. [19] and [21]. For 
^^Cr atoms, the natural DD/contact interaction strength 
ratio is less than 0.1 [l2[. Therefore, essential attenua- 
tion of the contact interaction by means of the FR is 
necessary [24]. One consequence of this requirement to 
the experiment is that three-body losses, which are also 
induced by the FR, of order 10"^^ cm^ s"^ [27], will set 
an upper limit on the free-evolution time and impose an 
intrinsic time dependence. Preliminary numerical sim- 
ulations including a quintic loss term of this magnitude 
indicate that soliton behavior is not considerably affected 
by these losses. The remaining difficulty (common with 
that for fundamental solitons predicted in Ref. p^]) is 
the necessity to invert the sign of the DD interaction by 
means of the rapidly rotating magnetic field fl7^ , a tech- 
nique which still has to be experimentally demonstrated. 

In conclusion, using the variational analysis and di- 
rect simulations of the GP equation in three dimensions, 
we have predicted the existence of stable quasi-2D vortex 
solitons with topological charge ^ = 1 in the dipolar BEG 
with atomic moments polarized perpendicular to the 2D 
plane, and the inverted sign of the dipole-dipole interac- 
tion. While energy minima exist also for solitary vortices 
with = 2, the resulting soliton is unstable and splits into 
a quadrupole set after the evolution in the course of a few 
trap periods. Future work will explore the possibility of 
anisotropic solitary vortices with an in-plane polarization 
axis and the natural sign of the DD interactions. 
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